/*
 * inverse_normal_distribution.cpp
 *
 *  Created on: Mar 13, 2009
 *      Author: genady
 */

#include <iostream>
#include <stdexcept>
#include <cmath>
#include <utility>
#include "inverse_normal_distribution.hpp"
#include "../utils/math_utils.hpp"
using namespace std;

static const double __a[] = {
	-3.969683028665376E+01,
	2.209460984245205E+02,
	-2.759285104469687E+02,
	1.383577518672690E+02,
	-3.066479806614716E+01,
	2.506628277459239E+00
};

static const double __b[] = {
	-5.447609879822406E+01,
	1.615858368580409E+02,
	-1.556989798598866E+02,
	6.680131188771972E+01,
	-1.328068155288572E+01
};

static const double __c[] = {
	-7.784894002430293E-03,
	-3.223964580411365E-01,
	-2.400758277161838E+00,
	-2.549732539343734E+00,
	4.374664141464968E+00,
	2.938163982698783E+00
};

static const double __d[] = {
	7.784695709041462E-03,
	3.224671290700398E-01,
	2.445134137142996E+00,
	3.754408661907416E+00
};

static const double __pLow = 0.02425;
static const double __pHigh = 1.0 - __pLow;

InverseNormalCDF::InverseNormalCDF( const double &mean, const double &rms ) throw()
	: _mean( mean ), _rms( rms ) {
}

InverseNormalCDF::InverseNormalCDF( const InverseNormalCDF &other ) throw()
	: _mean(other._mean), _rms(other._rms) {
}

InverseNormalCDF& InverseNormalCDF::operator= ( const InverseNormalCDF &other ) throw() {
	_mean = other._mean;
	_rms = other._rms;
	return *this;
}

bool InverseNormalCDF::operator== ( const InverseNormalCDF &other ) const throw() {
	return _mean == other._mean && _rms == other._rms;
}

double InverseNormalCDF::approximation4LowerRegion( const double &p ) const {
	double q = sqrt( -2*log( p ) );
	double x = (((((__c[0]*q + __c[1])*q + __c[2])*q + __c[3])*q + __c[4])*q + __c[5]) /
				((((__d[0]*q + __d[1])*q + __d[2])*q + __d[3])*q + 1);
	return x;
}

double InverseNormalCDF::approximation4CentralRegion( const double &p ) const {
	double q = p - 0.5;
	double r = q*q;
	double x = (((((__a[0]*r + __a[1])*r + __a[2])*r + __a[3])*r + __a[4])*r + __a[5])*q /
				(((((__b[0]*r + __b[1])*r + __b[2])*r + __b[3])*r + __b[4])*r + 1);
	return x;
}

double InverseNormalCDF::approximation4UpperRegion( const double &p ) const {
	double q = sqrt( -2*log( 1 - p ) );
	double x = -(((((__c[0]*q + __c[1])*q + __c[2])*q + __c[3])*q + __c[4])*q + __c[5]) /
				((((__d[0]*q + __d[1])*q + __d[2])*q + __d[3])*q + 1);
	return x;
}

double InverseNormalCDF::approximation( const double &p ) const throw(invalid_argument) {
	if( p > 0.0 && p < __pLow )
		return approximation4LowerRegion( p );
	else if( p >= __pLow && p <= __pHigh )
		return approximation4CentralRegion( p );
	else if( p >= __pHigh && p < 1.0 )
		return approximation4UpperRegion( p );
	else
		throw invalid_argument( "0 < confidence < 1" );
}

ConfidenceInterval::ConfidenceInterval( const double &mean, const double &rms, const double &confidence )
		throw(invalid_argument)
		: _mean(mean), _rms(rms), _confidence(confidence), _invNormCDF(mean, rms) {
	if( rms < 0 )
		throw invalid_argument( "rms >= 0" );
	double approximation = _invNormCDF.approximation( 0.5*(1.0 - confidence) );
	_lMargin = _mean + approximation*_rms;
	_rMargin = 2.0*_mean - _lMargin;
}

ConfidenceInterval::ConfidenceInterval( const ConfidenceInterval& other ) throw()
		: _mean(other._mean), _rms(other._rms), _confidence(other._confidence),
		_lMargin(other._lMargin), _rMargin(other._rMargin), _invNormCDF(other._invNormCDF)  {
}

ConfidenceInterval& ConfidenceInterval::operator= ( const ConfidenceInterval& other ) throw() {
	_mean = other._mean;
	_rms = other._rms;
	_confidence = other._confidence;
	_invNormCDF = other._invNormCDF;
	return *this;
}

bool ConfidenceInterval::operator== ( const ConfidenceInterval& other ) const throw() {
	return _invNormCDF == other._invNormCDF && _confidence == other._confidence;
}

pair<double, double> ConfidenceInterval::getInterval() const {
	return exp( getLeftMargin(), getRightMargin() );
}

ConfidenceInterval& ConfidenceInterval::recalculate( const double &confidence )
		throw(invalid_argument) {
	_confidence = confidence;
	double approximation = _invNormCDF.approximation( 0.5*(1.0 - confidence) );
	_lMargin = _mean + approximation*_rms;
	_rMargin = 2.0*_mean - _lMargin;
	return *this;
}

std::pair<double, double> ConfidenceInterval::getMargins() const
{
	return std::make_pair(_lMargin, _rMargin);
}

